
#### 4.1 MODEL OUTPUT LEGITIMATION TECHNOCRATIC ####

pacman::p_load(dplyr, readr, stringr, tidytext, quanteda, SnowballC, wordcloud2, purrr, fixest, tidyr, modelsummary, kableExtra)

source("./4_0_Models_output_data_generation.R")

papers <- readRDS("./data/scopus_results.rds") %>%
  select(dc_identifier, abstract) %>%
  unnest_sentences(input = abstract, output = abstract) %>%
  mutate(filterword = str_detect(abstract, "best practice|rational|scientific|technical|professional|neutral|knowledge|quality|innovation|solution|solve|efficient|efficiency|compatible|compatibility|effective|coordinate|coordination|streamline|strategic|systematic|optimize|optimal|optimization|cost-effective|cost-efficient|quality|performance|expert|experts|expertise|functional")) %>%
  filter(filterword == TRUE) %>%
  select(-filterword)

paper_topwords <- papers %>%
  select(text = abstract) %>%
  mutate(text = str_remove_all(text, "[0-9]+")) %>%
  mutate(text = str_remove_all(text, "[[:punct:]]")) %>%
  mutate(text = str_remove_all(text, symbols)) %>%
  unnest_tokens(token, text, token = "words") %>%
  anti_join(tidytext::stop_words %>% rename(token = word), by = "token") %>%
  mutate(token = wordStem(token, language = "english")) %>%
  count(token, sort = TRUE) %>%
  ungroup() %>%
  arrange(desc(n))

paper_words <- paper_topwords %>%
  distinct(token, .keep_all = TRUE) %>%
  filter(n > 500) 

wordcloud2(data = paper_words, size = 1, shape = 'pentagon')

tech_ling1 <- tibble(token = c(paper_words %>% pull(token)),
                      dimension = "papers")

# saveRDS(tech_ling1, file = "./data/technocratic_dict_papers.rds")

## PatentView

class_patents <- readRDS("./data/class_patents.rds")

patent_sentence <- class_patents %>%
  filter(patent_year > 2003 & patent_year < 2023) %>%
  distinct(patent_number, .keep_all = TRUE) %>%
  select(patent_number , patent_abstract) %>%
  unnest_sentences(input = patent_abstract, output = patent_abstract)

patents <- patent_sentence %>%
  mutate(filter_word = str_detect(patent_abstract, "scientific|technical|professional|neutral|knowledge|quality|innovation|solution|solve|efficient|efficiency|compatible|compatibility|effective|coordinate|coordination|streamline|strategic|systematic|optimize|optimal|optimization|cost-effective|cost-efficient|quality|performance|expert|experts|expertise|functional")) %>%
  filter(filter_word == TRUE) %>%
  select(-filter_word) %>%
  select(text = patent_abstract) %>%
  mutate(text = str_remove_all(text, "[0-9]+")) %>%
  mutate(text = str_remove_all(text, "[[:punct:]]")) %>%
  mutate(text = str_remove_all(text, symbols)) %>%
  unnest_tokens(token, text, token = "words") %>%
  anti_join(tidytext::stop_words %>% rename(token = word), by = "token") %>%
  mutate(token = wordStem(token, language = "english")) %>%
  count(token, sort = TRUE) %>%
  ungroup() %>%
  arrange(desc(n))

patents_words <- patents %>%
  distinct(token, .keep_all = TRUE) %>%
  filter(n > 15000) 

wordcloud2(data = patents_words, size = 1, shape = 'pentagon')

tech_ling2 <- tibble(token = c(patents_words %>% pull(token)),
                       dimension = "patent")

# saveRDS(tech_ling2, file = "./data/technocratic_dict_patents.rds")

## Merge with news text and generate shares

standard_news_token <- standard_news_text %>%
  unnest_tokens(token, text, token = "words") %>%
  anti_join(stopwords %>% rename(token = word), by = "token") 

paper_share <- standard_news_token %>%
  mutate(token = wordStem(token, language = "english")) %>%
  mutate(tech_ling1 = ifelse(token %in% c(tech_ling1$token), 1, 0)) %>%
  group_by(tech_ling1, rowid, year) %>%
  count() %>%
  spread(tech_ling1, n) %>%
  mutate(share_paper = `1`/(`0`+`1`)) %>%
  ungroup() %>%
  select(-`1`, -`0`) %>%
  left_join(news_text %>% select(rowid, name) %>% unique(), by = join_by(rowid))

patent_share <- standard_news_token %>%
  mutate(token = wordStem(token, language = "english")) %>%
  mutate(tech_ling2 = ifelse(token %in% c(tech_ling2$token), 1, 0)) %>%
  group_by(tech_ling2, rowid, year) %>%
  count() %>%
  spread(tech_ling2, n) %>%
  mutate(share_patent = `1`/(`0`+`1`)) %>%
  ungroup() %>%
  select(-`1`, -`0`) %>%
  left_join(news_text %>% select(rowid, name) %>% unique(), by = join_by(rowid))

standard_news_text_shares <- standard_news_text %>%
  left_join(paper_share, by = join_by(year, name, rowid)) %>%
  left_join(patent_share, by = join_by(year, name, rowid)) %>%
  unique()

mean_share_techleg <- standard_news_text_shares %>%
  group_by(name) %>%
  summarise(mean_share_paper = mean(share_paper, na.rm = TRUE),
            mean_share_patent = mean(share_patent, na.rm = TRUE))

output_tech_leg <- standard_news_text_shares %>%
  inner_join(mean_share_techleg, by = join_by(name)) %>%
  ungroup() %>%
  distinct(name, stdno, year, date, .keep_all = TRUE)

output_tech_leg %>%
  group_by(physical, societal) %>%
  count()

## Models

mod1 <- feols(mean_share_paper ~ physical | year + committee, 
              cluster = "year+committee",
              data = output_tech_leg) 

mod2 <- feols(mean_share_patent ~ physical | year + committee, 
              cluster = "year+committee",
              data = output_tech_leg)

mod3 <- feols(share_technocratic ~ physical | year + committee, 
              cluster = "year+committee", 
              data = legitimacy_coding)

## Models with societal standards

mod1_s <- feols(mean_share_paper ~ societal | year + committee, 
              cluster = "year+committee",
              data = output_tech_leg)

mod2_s <- feols(mean_share_patent ~ societal | year + committee, 
              cluster = "year+committee",
              data = output_tech_leg)

mod3_s <- feols(share_technocratic ~ societal | year + committee,
                cluster = "year+committee",
                data = legitimacy_coding)

cm <- c("physical" = "Physical",
        "societal" = "Societal")

models <- list(mod1, mod2, mod3,
               mod1_s, mod2_s, mod3_s)

names(models) <- c("Research papers", "Patents", "Legitimation Statements",
                   "Research papers", "Patents", "Legitimation Statements")

mod <- modelsummary(models,
                    stars = TRUE,
                    gof_omit = 'AIC|BIC|Within|Std.Errors|RMSE|FE',
                    fmt = 3,
                    coef_map = cm,
                    #output = "latex",
                    statistic = c("{std.error}"),
                    notes = c("Fixed effects: Year and committee",
                              "Clustered standard errors by committee and year",
                              "Coverage: 1999 - 2022")) %>%
  kable_styling(font_size = 11, full_width = FALSE) %>%
  add_header_above(c(" " = 1, "Percent technocratic legitimation words/statements" = 6))

mod

## Coefficient plot

tech_leg1_physical <- broom::tidy(feols(mean_share_paper ~ physical | year + committee, 
                               cluster = "year+committee",
                             data = output_tech_leg %>%
                               mutate(mean_share_paper = scale(mean_share_paper, center = TRUE, scale = TRUE),
                                      physical = scale(physical, center = TRUE, scale = TRUE))),
                       conf.int = TRUE) %>% filter(term == "physical") %>% mutate(term = ifelse(term == "physical", "Research papers", NA))

tech_leg2_physical <- broom::tidy(feols(mean_share_patent ~ physical | year+committee, 
                               cluster = "year+committee",
                              data = output_tech_leg %>%
                                mutate(mean_share_patent = scale(mean_share_patent, center = TRUE, scale = TRUE),
                                       physical = scale(physical, center = TRUE, scale = TRUE))),
                        conf.int = TRUE) %>% filter(term == "physical") %>% mutate(term = ifelse(term == "physical", "Patents", NA))

tech_leg3_physical <- broom::tidy(feols(share_technocratic ~ physical | year+committee, 
                                        cluster = "year+committee",
                                        data = legitimacy_coding %>%
                                          mutate(share_technocratic = scale(share_technocratic, center = TRUE, scale = TRUE),
                                                 physical = scale(physical, center = TRUE, scale = TRUE))),
                                  conf.int = TRUE) %>% filter(term == "physical") %>% mutate(term = ifelse(term == "physical", "Legitimation  statements", NA))

tech_leg1_societal <- broom::tidy(feols(mean_share_paper ~ societal | year + committee, 
                               cluster = "year+committee",
                               data = output_tech_leg %>%
                                 mutate(mean_share_paper = scale(mean_share_paper, center = TRUE, scale = TRUE),
                                        societal = scale(societal, center = TRUE, scale = TRUE))),
                         conf.int = TRUE) %>% filter(term == "societal") %>% mutate(term = ifelse(term == "societal", "Research papers", NA))

tech_leg2_societal <- broom::tidy(feols(mean_share_patent ~ societal | year+committee, 
                               cluster = "year+committee",
                               data = output_tech_leg %>%
                                 mutate(mean_share_patent = scale(mean_share_patent, center = TRUE, scale = TRUE),
                                        societal = scale(societal, center = TRUE, scale = TRUE))),
                         conf.int = TRUE) %>% filter(term == "societal") %>% mutate(term = ifelse(term == "societal", "Patents", NA))

tech_leg3_societal <- broom::tidy(feols(share_technocratic ~ societal | year+committee, 
                                        cluster = "year+committee",
                                        data = legitimacy_coding %>%
                                          mutate(share_technocratic = scale(share_technocratic, center = TRUE, scale = TRUE),
                                                 societal = scale(societal, center = TRUE, scale = TRUE))),
                                  conf.int = TRUE) %>% filter(term == "societal") %>% mutate(term = ifelse(term == "societal", "Legitimation  statements", NA))

coef_plot_output_tech <- bind_rows(tech_leg1_physical %>% mutate(dv = "Research papers", model = "Technocratic output legitimation", standard = "physical"),
                                   tech_leg2_physical %>% mutate(dv = "Patents", model = "Technocratic output legitimation", standard = "physical"),
                                   tech_leg3_physical %>% mutate(dv = "Legitimation  statements", model = "Technocratic output legitimation", standard = "physical")) %>%
  bind_rows(tech_leg1_societal %>% mutate(dv = "Research papers", model = "Technocratic output legitimation", standard = "societal"),
            tech_leg2_societal %>% mutate(dv = "Patents", model = "Technocratic output legitimation", standard = "societal"),
            tech_leg3_societal %>% mutate(dv = "Legitimation  statements", model = "Technocratic output legitimation", standard = "societal"))


#### ROBUSTNESS ####

# ### 200 / 10000 ###
# 
# paper_words <- paper_topwords %>%
#   distinct(token, .keep_all = TRUE) %>%
#   filter(n > 200) 
# 
# patents_words <- patents %>%
#   distinct(token, .keep_all = TRUE) %>%
#   filter(n > 10000)
# 
# tech_ling1 <- tibble(token = c(paper_words %>% pull(token)),
#                      dimension = "papers")
# 
# tech_ling2 <- tibble(token = c(patents_words %>% pull(token)),
#                      dimension = "patent")
# 
# standard_news_token <- standard_news_text %>%
#   unnest_tokens(token, text, token = "words") %>%
#   anti_join(stopwords %>% rename(token = word), by = "token") 
# 
# paper_share <- standard_news_token %>%
#   mutate(token = wordStem(token, language = "english")) %>%
#   mutate(tech_ling1 = ifelse(token %in% c(tech_ling1$token), 1, 0)) %>%
#   group_by(tech_ling1, rowid, year) %>%
#   count() %>%
#   spread(tech_ling1, n) %>%
#   mutate(share_paper = `1`/(`0`+`1`)) %>%
#   ungroup() %>%
#   select(-`1`, -`0`) %>%
#   left_join(news_text %>% select(rowid, name) %>% unique(), by = join_by(rowid))
# 
# patent_share <- standard_news_token %>%
#   mutate(token = wordStem(token, language = "english")) %>%
#   mutate(tech_ling2 = ifelse(token %in% c(tech_ling2$token), 1, 0)) %>%
#   group_by(tech_ling2, rowid, year) %>%
#   count() %>%
#   spread(tech_ling2, n) %>%
#   mutate(share_patent = `1`/(`0`+`1`)) %>%
#   ungroup() %>%
#   select(-`1`, -`0`) %>%
#   left_join(news_text %>% select(rowid, name) %>% unique(), by = join_by(rowid))
# 
# standard_news_text_shares <- standard_news_text %>%
#   left_join(paper_share, by = join_by(year, name, rowid)) %>%
#   left_join(patent_share, by = join_by(year, name, rowid)) %>%
#   unique()
# 
# mean_share_techleg <- standard_news_text_shares %>%
#   group_by(name) %>%
#   summarise(mean_share_paper = mean(share_paper, na.rm = TRUE),
#             mean_share_patent = mean(share_patent, na.rm = TRUE))
# 
# output_tech_leg <- standard_news_text_shares %>%
#   inner_join(mean_share_techleg, by = join_by(name)) %>%
#   ungroup() %>%
#   distinct(name, stdno, year, date, .keep_all = TRUE)
# 
# mod1 <- feols(mean_share_paper ~ physical | year + committee, 
#               cluster = "year+committee",
#               data = output_tech_leg) 
# 
# mod2 <- feols(mean_share_patent ~ physical | year + committee, 
#               cluster = "year+committee",
#               data = output_tech_leg)
# 
# cm <- c("physical" = "Physical",
#         "societal" = "Societal")
# 
# models <- list(mod1, mod2)
# 
# names(models) <- c("Research papers", "Patents")
# 
# mod_rob1 <- modelsummary(models,
#                     stars = TRUE,
#                     gof_omit = 'AIC|BIC|Within|Std.Errors|RMSE|FE',
#                     fmt = 3,
#                     coef_map = cm,
#                     output = "latex",
#                     statistic = c("{std.error}"),
#                     notes = c("Fixed effects: Year and committee",
#                               "Clustered standard errors by committee and year",
#                               "Coverage: 1999 - 2022")) %>%
#   kable_styling(font_size = 11, full_width = FALSE) %>%
#   add_header_above(c(" " = 1, "Cutoff: 200" = 1, "Cutoff: 10000" = 1)) %>%
#   add_header_above(c(" " = 1, "Percent technocratic legitimation words" = 2))
# 
# ### 400 / 12000 ###
# 
# paper_words <- paper_topwords %>%
#   distinct(token, .keep_all = TRUE) %>%
#   filter(n > 400) 
# 
# patents_words <- patents %>%
#   distinct(token, .keep_all = TRUE) %>%
#   filter(n > 12000)
# 
# tech_ling1 <- tibble(token = c(paper_words %>% pull(token)),
#                      dimension = "papers")
# 
# tech_ling2 <- tibble(token = c(patents_words %>% pull(token)),
#                      dimension = "patent")
# 
# standard_news_token <- standard_news_text %>%
#   unnest_tokens(token, text, token = "words") %>%
#   anti_join(stopwords %>% rename(token = word), by = "token") 
# 
# paper_share <- standard_news_token %>%
#   mutate(token = wordStem(token, language = "english")) %>%
#   mutate(tech_ling1 = ifelse(token %in% c(tech_ling1$token), 1, 0)) %>%
#   group_by(tech_ling1, rowid, year) %>%
#   count() %>%
#   spread(tech_ling1, n) %>%
#   mutate(share_paper = `1`/(`0`+`1`)) %>%
#   ungroup() %>%
#   select(-`1`, -`0`) %>%
#   left_join(news_text %>% select(rowid, name) %>% unique(), by = join_by(rowid))
# 
# patent_share <- standard_news_token %>%
#   mutate(token = wordStem(token, language = "english")) %>%
#   mutate(tech_ling2 = ifelse(token %in% c(tech_ling2$token), 1, 0)) %>%
#   group_by(tech_ling2, rowid, year) %>%
#   count() %>%
#   spread(tech_ling2, n) %>%
#   mutate(share_patent = `1`/(`0`+`1`)) %>%
#   ungroup() %>%
#   select(-`1`, -`0`) %>%
#   left_join(news_text %>% select(rowid, name) %>% unique(), by = join_by(rowid))
# 
# standard_news_text_shares <- standard_news_text %>%
#   left_join(paper_share, by = join_by(year, name, rowid)) %>%
#   left_join(patent_share, by = join_by(year, name, rowid)) %>%
#   unique()
# 
# mean_share_techleg <- standard_news_text_shares %>%
#   group_by(name) %>%
#   summarise(mean_share_paper = mean(share_paper, na.rm = TRUE),
#             mean_share_patent = mean(share_patent, na.rm = TRUE))
# 
# output_tech_leg <- standard_news_text_shares %>%
#   inner_join(mean_share_techleg, by = join_by(name)) %>%
#   ungroup() %>%
#   distinct(name, stdno, year, date, .keep_all = TRUE)
# 
# mod1 <- feols(mean_share_paper ~ physical | year + committee, 
#               cluster = "year+committee",
#               data = output_tech_leg) 
# 
# mod2 <- feols(mean_share_patent ~ physical | year + committee, 
#               cluster = "year+committee",
#               data = output_tech_leg)
# 
# cm <- c("physical" = "Physical",
#         "societal" = "Societal")
# 
# models <- list(mod1, mod2)
# 
# names(models) <- c("Research papers", "Patents")
# 
# mod_rob2 <- modelsummary(models,
#                          stars = TRUE,
#                          gof_omit = 'AIC|BIC|Within|Std.Errors|RMSE|FE',
#                          fmt = 3,
#                          coef_map = cm,
#                          output = "latex",
#                          statistic = c("{std.error}"),
#                          notes = c("Fixed effects: Year and committee",
#                                    "Clustered standard errors by committee and year",
#                                    "Coverage: 1999 - 2022")) %>%
#   kable_styling(font_size = 11, full_width = FALSE) %>%
#   add_header_above(c(" " = 1, "Cutoff: 400" = 1, "Cutoff: 12000" = 1)) %>%
#   add_header_above(c(" " = 1, "Percent technocratic legitimation words" = 2))
# 
# ### 800 / 20000 ###
# 
# paper_words <- paper_topwords %>%
#   distinct(token, .keep_all = TRUE) %>%
#   filter(n > 800) 
# 
# patents_words <- patents %>%
#   distinct(token, .keep_all = TRUE) %>%
#   filter(n > 20000)
# 
# tech_ling1 <- tibble(token = c(paper_words %>% pull(token)),
#                      dimension = "papers")
# 
# tech_ling2 <- tibble(token = c(patents_words %>% pull(token)),
#                      dimension = "patent")
# 
# standard_news_token <- standard_news_text %>%
#   unnest_tokens(token, text, token = "words") %>%
#   anti_join(stopwords %>% rename(token = word), by = "token") 
# 
# paper_share <- standard_news_token %>%
#   mutate(token = wordStem(token, language = "english")) %>%
#   mutate(tech_ling1 = ifelse(token %in% c(tech_ling1$token), 1, 0)) %>%
#   group_by(tech_ling1, rowid, year) %>%
#   count() %>%
#   spread(tech_ling1, n) %>%
#   mutate(share_paper = `1`/(`0`+`1`)) %>%
#   ungroup() %>%
#   select(-`1`, -`0`) %>%
#   left_join(news_text %>% select(rowid, name) %>% unique(), by = join_by(rowid))
# 
# patent_share <- standard_news_token %>%
#   mutate(token = wordStem(token, language = "english")) %>%
#   mutate(tech_ling2 = ifelse(token %in% c(tech_ling2$token), 1, 0)) %>%
#   group_by(tech_ling2, rowid, year) %>%
#   count() %>%
#   spread(tech_ling2, n) %>%
#   mutate(share_patent = `1`/(`0`+`1`)) %>%
#   ungroup() %>%
#   select(-`1`, -`0`) %>%
#   left_join(news_text %>% select(rowid, name) %>% unique(), by = join_by(rowid))
# 
# standard_news_text_shares <- standard_news_text %>%
#   left_join(paper_share, by = join_by(year, name, rowid)) %>%
#   left_join(patent_share, by = join_by(year, name, rowid)) %>%
#   unique()
# 
# mean_share_techleg <- standard_news_text_shares %>%
#   group_by(name) %>%
#   summarise(mean_share_paper = mean(share_paper, na.rm = TRUE),
#             mean_share_patent = mean(share_patent, na.rm = TRUE))
# 
# output_tech_leg <- standard_news_text_shares %>%
#   inner_join(mean_share_techleg, by = join_by(name)) %>%
#   ungroup() %>%
#   distinct(name, stdno, year, date, .keep_all = TRUE)
# 
# mod1 <- feols(mean_share_paper ~ physical | year + committee, 
#               cluster = "year+committee",
#               data = output_tech_leg) 
# 
# mod2 <- feols(mean_share_patent ~ physical | year + committee, 
#               cluster = "year+committee",
#               data = output_tech_leg)
# 
# cm <- c("physical" = "Physical",
#         "societal" = "Societal")
# 
# models <- list(mod1, mod2)
# 
# names(models) <- c("Research papers", "Patents")
# 
# mod_rob3 <- modelsummary(models,
#                          stars = TRUE,
#                          gof_omit = 'AIC|BIC|Within|Std.Errors|RMSE|FE',
#                          fmt = 3,
#                          coef_map = cm,
#                          output = "latex",
#                          statistic = c("{std.error}"),
#                          notes = c("Fixed effects: Year and committee",
#                                    "Clustered standard errors by committee and year",
#                                    "Coverage: 1999 - 2022")) %>%
#   kable_styling(font_size = 11, full_width = FALSE) %>%
#   add_header_above(c(" " = 1, "Cutoff: 800" = 1, "Cutoff: 20000" = 1)) %>%
#   add_header_above(c(" " = 1, "Percent technocratic legitimation words" = 2))
# 
# ### 1200 / 25000 ###
# 
# paper_words <- paper_topwords %>%
#   distinct(token, .keep_all = TRUE) %>%
#   filter(n > 1200) 
# 
# patents_words <- patents %>%
#   distinct(token, .keep_all = TRUE) %>%
#   filter(n > 25000)
# 
# tech_ling1 <- tibble(token = c(paper_words %>% pull(token)),
#                      dimension = "papers")
# 
# tech_ling2 <- tibble(token = c(patents_words %>% pull(token)),
#                      dimension = "patent")
# 
# standard_news_token <- standard_news_text %>%
#   unnest_tokens(token, text, token = "words") %>%
#   anti_join(stopwords %>% rename(token = word), by = "token") 
# 
# paper_share <- standard_news_token %>%
#   mutate(token = wordStem(token, language = "english")) %>%
#   mutate(tech_ling1 = ifelse(token %in% c(tech_ling1$token), 1, 0)) %>%
#   group_by(tech_ling1, rowid, year) %>%
#   count() %>%
#   spread(tech_ling1, n) %>%
#   mutate(share_paper = `1`/(`0`+`1`)) %>%
#   ungroup() %>%
#   select(-`1`, -`0`) %>%
#   left_join(news_text %>% select(rowid, name) %>% unique(), by = join_by(rowid))
# 
# patent_share <- standard_news_token %>%
#   mutate(token = wordStem(token, language = "english")) %>%
#   mutate(tech_ling2 = ifelse(token %in% c(tech_ling2$token), 1, 0)) %>%
#   group_by(tech_ling2, rowid, year) %>%
#   count() %>%
#   spread(tech_ling2, n) %>%
#   mutate(share_patent = `1`/(`0`+`1`)) %>%
#   ungroup() %>%
#   select(-`1`, -`0`) %>%
#   left_join(news_text %>% select(rowid, name) %>% unique(), by = join_by(rowid))
# 
# standard_news_text_shares <- standard_news_text %>%
#   left_join(paper_share, by = join_by(year, name, rowid)) %>%
#   left_join(patent_share, by = join_by(year, name, rowid)) %>%
#   unique()
# 
# mean_share_techleg <- standard_news_text_shares %>%
#   group_by(name) %>%
#   summarise(mean_share_paper = mean(share_paper, na.rm = TRUE),
#             mean_share_patent = mean(share_patent, na.rm = TRUE))
# 
# output_tech_leg <- standard_news_text_shares %>%
#   inner_join(mean_share_techleg, by = join_by(name)) %>%
#   ungroup() %>%
#   distinct(name, stdno, year, date, .keep_all = TRUE)
# 
# mod1 <- feols(mean_share_paper ~ physical | year + committee, 
#               cluster = "year+committee",
#               data = output_tech_leg) 
# 
# mod2 <- feols(mean_share_patent ~ physical | year + committee, 
#               cluster = "year+committee",
#               data = output_tech_leg)
# 
# cm <- c("physical" = "Physical",
#         "societal" = "Societal")
# 
# models <- list(mod1, mod2)
# 
# names(models) <- c("Research papers", "Patents")
# 
# mod_rob4 <- modelsummary(models,
#                          stars = TRUE,
#                          gof_omit = 'AIC|BIC|Within|Std.Errors|RMSE|FE',
#                          fmt = 3,
#                          coef_map = cm,
#                          output = "latex",
#                          statistic = c("{std.error}"),
#                          notes = c("Fixed effects: Year and committee",
#                                    "Clustered standard errors by committee and year",
#                                    "Coverage: 1999 - 2022")) %>%
#   kable_styling(font_size = 11, full_width = FALSE) %>%
#   add_header_above(c(" " = 1, "Cutoff: 1200" = 1, "Cutoff: 25000" = 1)) %>%
#   add_header_above(c(" " = 1, "Percent technocratic legitimation words" = 2))
# 
# 
# mod_rob1
# mod_rob2
# mod_rob3
# mod_rob4

## Clean

rm(list=ls()[! ls() %in% c("coef_plot_input_dem", "coef_plot_input_tech", "coef_plot_output_dem", "coef_plot_output_tech",
                           "liaison_dat_merged_sh", "liaison_dat_merged_frac", "participation_dat_merged_sh", "participation_dat_merged_frac",
                           "tech_leg",
                           "output_tech_leg", "output_dem_leg",
                           "legitimacy_coding")])
gc()


